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We perform a numerical study of the critical regime at the threshold of black hole formation in 
the spherically symmetric, general relativistic collapse of collisionless matter. The coupled Einstein- 
Vlasov equations are solved using a particle-mesh method in which the evolution of the phase-space 
distribution function is approximated by a set of particles (or, more precisely, infinitesimally thin 
shells) moving along geodesies of the spacetime. Individual particles may have non-zero angular 
momenta, but spherical symmetry dictates that the total angular momentum of the matter distri- 
^ I ^ ' bution vanish. In accord with previous work by Rein et al ||l|], our results indicate that the critical 

^^ , behavior in this model is Type I; that is, the smallest black hole in each parametrized family has 

^N ■ a finite mass. We present evidence that the critical solutions are characterized by unstable, static 

spacetimes, with non-trivial distributions of radial momenta for the particles. As expected for Type 
I solutions, we also find power-law scaling relations for the lifetimes of near-critical configurations 
as a function of parameter-space distance from criticality. 
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I. INTRODUCTION 



^^ I Critical phenomena at the threshold of black hole formation were originally discovered in studies of the spherically 

C_^ ■ symmetric, general relativistic collapse of a minimally coupled scalar field [EJ. Similar behavior has now been found 
j^ in many different scenarios, including the collapse of gravitational waves, perfect fluids, Yang-Mills fields and scalar 
^— V . fields in anti-de Sitter spacetime (for a review see Q). Relatively little work has been done on the critical collapse 
of collisionless matter. To date, the only detailed study of the black-hole threshold in the Einstein- Vlasov model is 
due to Rein, Rendall and Schaeffer pi. These authors found evidence that for spherically symmetric collapse with 
non-zero angular momenta distributions, the threshold black hole mass is finite (Type I behavior). In this paper we 
^H, summarize the results of M which corroborate and extend the previous work of Rein et al. 
I ■ The paper is organized as follows: In section |l|, we outline the specific form of the Einstein- Vlasov equations we 

Kr\' have solved, and make some contact with the particle-mesh (PM) method which is subsequently used to numerically 
solve these equations. Here we follow the approach of Shapiro and Tcukolsky g, which has been successfully used 

describes our 



to model the dynamics of spherically-symmetric, relativistic clusters of stars If, m, M. Section III 



numerical techniques per se, while section IV contains our main results, including evidence that the critical solutions 



\^ in this model are characterized by static geometries and satisfy the type of scaling expected of Type I solutions. 



Finally, some brief concluding remarks are made in section M. 

We use geometric units, G = c = 1, throughout the paper. Abstract spacetime indices are generally denoted by 
o and b, while /it, v and fc, I are used for spacetime and spatial component indices, respectively. Finally, subscript j's 
label specific particles, while subscript j's are generally used for finite-difference indexing. 

II. FORMALISM AND EQUATIONS OF MOTION 

The dynamical state of collisionless matter can be described by a distribution function, f{x°',pa)'. 

![x\pa) = dN/dVj, (I) 
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where N is the particle number and Vp is the phase-space volume. In the current case, the volume in phase space is 
conserved during the evolution of the system (Liouville's theorem). This implies that the distribution function is also 
a conserved quantity: 

df{t,x'',pk) _ , . 

dt " " ^ > 

This is the collisionless Boltzmann, or Vlasov, equation. This equation, coupled to Einstein's equations. Gab = ^T^Tab, 
all restricted to spherical symmetry, i.e.: 

fit,x'',pk) = fit,Rx\Rpk) with Re 30(3) fc = l,2,3 (3) 

form the system we wish to solve numerically. 

A. Maximal-areal coordinate system. 

As with any problem in "3+1" ("space + time") numerical relativity, we want to specify initial data on a spatial 
hypersurface and then evolve these data in time. To do this, we need to split the Einstein equations into a set of 
constraint equations (equations that must be satisfied at each instant of time) and dynamical or evolution equations 
(equations that tell us how to evolve the geometric quantities in time). We carry out this splitting using the 3+1 
formalism due to Arnowitt, Deser and Misner (ADM) (for reviews of this formalism, see g and flCJ]). 

We restrict attention to spherical symmetry and adopt coordinates (t, r, 6, ip) with the usual spherical-polar topology. 
We are left with the freedom to choose our radial and time coordinates, and have chosen maximal-areal coordinates. 
As the name suggests, in this system the radial coordinate is areal, so that the proper area of 2-spheres with radius r 
is 47rr^. The time coordinate is fixed by demanding that the t = constant, 3-slices be maximal, i.e. that the trace of 
the extrinsic curvature, K{t,r) = K^i{t,r), identically vanish on each slice. This leads to a slicing condition on the 
lapse Junction, a {t, r), which must be satisfied at each instant of time. 

With these choices, the spacetime metric takes the specific form: 

ds^ = {-a{t, rf + a{t, rf[3{t, rf) dt^ + 2a^/3 dtdr + a^ dr^ + r^ {dO^ + sin^ Odip^) (4) 

where P{t,r) is the radial component of the shift vector, (3^ — (/?, 0, 0). A sufhcient set of equations for determining 
the geometric quantities, a{t,r), K^g[t,r), a{t,r) and f3{t,r) is then: 
Hamiltonian Constraint: 

— ^a'^rK%^ + 'i-Kra'^p+ — il~a'^) (5) 

a 2 2r 



Momentum, Constraint: 



Slicing Condition: 



Areal Coordinate Condition: 



K%' = —K%~Anjr (6) 
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^ ( 2r- + a^ - 1 ) + A^a'a {S - 3p) (7) 



/3 = arK\ (8) 



Here, ' denotes differentiation with respect to r, and the last formula is derived from dtK — 0, using K = Q. In 
addition, p(t, r), jr{t, r) and S{t, r), which are discussed in detail in the next section, are the local energy density, the 
local current density and the trace of the spatial part of the stress-energy density, respectively. We note that we have 
chosen to implement a fully constrained evolution, which in this case means that we use the constraint equations, 
rather than evolution equations, to update a and K^e- 

During our simulations, we also compute the mass aspect function, M{t,r): 

2 I a2 ^2 2 



M{t, r) = ^(l + ^-4UMl + r'K\' - ^ ) , (9) 



which, among other useful diagnostic purposes, allows us detect the formation of apparent horizons. Specifically, 
when 2M{t,r)/r = 1 — l/a^ + 0^ jo? , becomes equal to 1, a marginally trapped surface has been formed. We can see 
this by computing the expansion of the outgoing null geodesies (see for instance [0), which in these coordinates can 
be written as: 

l-a(t,r)ri^Vt,r)-l-^^^^. (10) 

a(t,r) 

Therefore, if the outgoing expansion is zero, Xja} — 0^ /a^, and 2M{t,r)/r — 1. 

B. Stress-energy tensor 

In this section we explain how we calculate the stress-energy quantities p(t,r), jr{t,r) and S(t,r) that appear 
in equations (&p|). Adopting a Monte Carlo approach, we approximate the distribution function (1^) by a set of N 
"spherical particles" , which actually represent infinitesimally thin spherical shells of matter. Since these particles only 
interact with each other gravitationally, we have 

N 

T'"' — \^ Xf^'^ (11) 

4 = 1 

where Tf^'^ is the stress energy tensor for a single particle. For a point particle we have 

Tr^^^S{r-n{t)), (12) 

where, pf are the components of the 4- momentum of the i-th particle, rrii is its rest mass, ri{t) is its radial position 
at time t, and S is the usual Dirac (5-function. In maximal-areal coordinates, the single-particle contributions to the 
quantities p, S and jr then take the form: 

[Ph =«'[r"], (13) 

1 .„ , 1 

[jr\^a[T\]^. (15) 



[S\ = - [Trrl. + 3 [Tee], + ^-^^ [T,^], (14) 



We now relax the point particle approximation and assume that each particle is a spherically symmetric shell of mass, 
uniformly distributed over a region Ar in radius (subsequently, Ar will be identified with the mesh spacing, h, used 
in the finite-difference solution of the geometrical equations). Each shell of matter is to be interpreted as an average 
over an ensemble of shells, each centered at r = r^, and with angular momentum vectors which point in all possible 
directions. Thus, for any shell, the net angular momentum is zero, I = 0, but |f|'^ = ^^ 7^ 0. The proper volume 
occupied by each particle is then given by: 

t r 2 i 

V, = —Ary^d(pde^A'n:Araa^-!-^, (16) 

rui J rrii 

and we can approximate the delta function that appears in dl2) by 1/Vi. This yields: 

(17) 
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Here [p*] ■ is defined by [p*]^ = a [p*]^, and [l]^ = [pg],. + [p^]^ /sin^ 6i is the square of the magnitude of the angular 



momentum of the z-th particle. The geometric quantities a and a are evaluated at r = r^ as described in Sec. Ill C 
Note that we have also defined 









(20) 



i.e., for S"" a the index a is summed over the angular coordinates. We then introduce quantities which do not explicitly 



depend on the geometrical quantities: [p\^ = a[p]j, \S^r 



[S\\, [S\]^ = a[SW and [>], = a[jr\. In our 



numerical implementation of the equations of motion, these definitions provide a clean separation of the particle 
updates and the updates of the geometry variables. 

We interpolate the one-particle quantities to the continuum and sum over all the particles to find the total values: 



N 



/» = E/'^(^ -''*)' 



(21) 



4=1 



where fi is any of the single-particle barred quantities defined above , / is the corresponding continuum quantity, and 
W{r — Ti) is an interpolation function defined in detail in Sec. IIIC (see equation (pi])). Having defined ( |2l] ) we can 
now write equations (|[^) as 
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(22) 

(23) 

(24) 
(25) 



C. Evolution equations 

Because there are no explicit interactions between the particles, their equations of motion are just the spacetime 
geodesic equations (the characteristics of the Vlasov equation) . These can be derived from the formula for parallel 
transport of a particle's four-momentum along its world line: 



It proves useful to recast these equations in terms of the quantities, Pr, p* and F 
p^ in terms of these variables as: 

p^{t,r)^P^-P{t,r)p\t,r) 
a'^lt, r) 



(26) 
Pe^ +P>fi'^/ sin^ 9. We can express 



(27) 



To compute total derivatives with respect to coordinate time we use 

d d dr dr d d p^ d 
dt dt dr dt dr dt p* 9r' 



(28) 



where here, and in the remainder of this section, r is the particle's proper time. Applying this operator to equation 
(|^) we get 



dp'' 
~dt 



1 dpr p^ f da p^ da 



a2 dt 



dt p* dr 



dp* ,fd(i p'' dp 



dt 



dt p* dr J 



(29) 



Substituting equations (27) and (G9) into equation (Gq) we obtain: 

da 



dpr 
~dt 



'a—-p 
dr 



dp 
dr 



Pr 



1 dapr I 

q3 Qj. pt pt j,3 ' 



(30) 



which is the evolution equation for pr- To derive the evolution equation for r, we use the definition of p'' (jj^ = dr/dr), 
which after some manipulation yields: 



The time component of the 4- momentum, p*, is calculated using the normalization condition p^Pp, — —m?: 



ap' = Jm^ + ^ + -^ (32) 

It is also convenient, as previously mentioned, to use p* = ap* rather than p* itself. Using this definition in equa- 
tions (pQ), ( pi] ) and d32) yields the final form of the particle equations of motion: 

dpr da_^ dp adapr'^ I'^a 

at or or a-^ or p^ p^r-^ 

dr ap> 



dt a^p^ 



(3 (34) 



p-* = \/m2 + ^ + ii. (35) 



III. NUMERICAL APPROACH 



As discussed in the previous section, we have adopted a Monte Carlo, particle-based strategy to the solution of the 
Vlasov equation. In this approach we generate an A^-particle sample of some specified initial distribution function 
f{0,x'',Pk), and then use dynamical evolution of the N particles to approximate the full dynamics of f{t,x'',pk)- The 
continuum limit is recovered in the limit N ^ oo and, in the absence of any sophisticated "importance sampling" 
techniques, we expect the level of statistical error in our particle calculations to be of the order of 1/yN. We couple 
the particles to the gravitational field by introducing a finite-difference mesh on which we approximately solve the 
geometric equations, and by introducing transfer operators which allow us to produce mesh-based representations of 
particle quantities and vice versa. "Particle-Mesh" , or PM, methods such as ours are commonly used in the solution 
of Boltzmann equations, particularly those involving long-range interactions, and the reader is referred to p^ for a 
detailed review of such techniques. 

Here we simply note that a PM method is generically characterized by the splitting of each discrete time step, 
t" -^ f"+^ into two stages: 1) the solution of the field equations on a finite- difference mesh, and 2) the updating of 
particle positions via discrete versions of their equations of motion. In our case, and as described in the previous 
section, at each time step the stress-energy quantities are calculated by considering each particle to be "smoothed" 
over a finite volume. 

A. The field equations 

We first explain how the equations ( p2[]25| ) for the geometry are solved numerically, assumi ng tha t we know the 
quantities p, jr, S''r, 6'°a (our computation of the stress-energy quantities is described in Section niC| ). The first two 



equations, equations (£2123), are integrated from the origin, r = 0, using the Isoda |13| integrator. The boundary 
conditions are given by the spherical symmetry of the spacetime, and by the demand that the spacetime be locally 
fiat at r = 0. They are a(i,0) = 1, and K\{t,0) = 0. 

We compute the values of the functions aj and K^ej on a uniform grid of N^. points, rj = (j — l)h, j = 1, ■ ■ ■ N^, 
where h = Ar — r^iax/iNr — 1), and r — r^ax is the outer edge of the computational domain. 

In order to compute the values at r = rj+i, we supply to Isoda the values of the functions at r = rj and the 
derivatives computed using equations (|22[]23|) at r = rj+1/2, using the average of p and jr at rj and r^+i. 

[ph+1/2 = \ {\p\j + \p\,+i) (36) 

[>],+i/2 - \ ([>], + [>],+i) . (37) 

Once we have calculated a, we can solve the slicing equation: 



a" = a'(^-'^ + '^(a'-\ + 2r^ + 4^aa {S\ + S^ - ip) , (38) 

with the boundary conditions: 

a'(i,0) = 0, (39) 

a{t,oo) = \. (40) 

Here the first condition follows from the demand that the slicing be regular at r = 0, and the second one follows 
from asymptotic flatness, plus the demand that t measure proper time at infinity. We solve (pq) using a second-order 
finite-difference approximation on the finite-difference mesh: 



h"^ 2h \ 2haj tj 

^ [ a,2 - 1 + 2r. ^+' ^ 



,2 \-= _ --^ 2ha, 
AT:a,a, ( i-^ + [Sf] ^ - 3 [pV ) (41) 



Jj 



Rearranging this equation gives us: 



1 , /A _ f^.A^.f^ fj 



where: 



/? + 2t"'-'-b?+*"'+F- a "'-' = '' <''2) 



Z'^^^^ls?^-! («' 



2 („2 1 , o„ «J+i-«J-A , .„„Y['^'''']j 



''J 



^. = fl «/ - 1 + ^^ Po,:' + 4™, L-^ + [5:]^. - 3 [p]^ (44) 



2/iaj 
In addition to (E2[) we have the boundary equation at ?' = 






which can be derived from the 0{h^) forward finite-difference approximation to a' = at r = 

— 3ai 4- 4a2 — a^ 



2h 
and equation (E3) with j = 2. We have also the boundary condition at r = r 



= , (46) 



where Af is the mass aspect function defined by equation (H). This approximation follows from the known represen- 



tation of the asymptotically-fiat Schwarzschild solution in maximal-areal coordinates. Equations (42), (45) and ( |47|) 
constitute a linear tridiagonal system that can be solved using a tridiagonal solver (we have used the LAPACK p4| ] 
routine dgtsv). 



B. The evolution equations 

To evolve the particles' positions and momenta we integrate the geodesic equations (p^p^. The values of the 
coefficients in these equations (basically products and quotients of a, a, (3, a', a' and /?') must be calculated at the 



particle positions, r^, using the values obtained at the mesh points, Tj. The mesh values are interpolated to the 
particle positions using the same operator kernel used to produce mesh values from particle quantities (this procedure 
is explained in the next section). The geodesic equations are also integrated using the Isoda routine. At discrete time 
t = t", given a particle's position, r", and radial momentum, p", we calculate the new position, r""*"^, and momentum, 
p"+^, at i = t"^^ = t" + At by supplying to Isoda the values of the metric functions and their spatial derivatives 
evaluated at t = t". Because we use the t — t" values of the geometric quantities in the particle updates, rather than, 
for example, values at i = i"+^/^, we expect our solution of the particle equations to have accuracy 0(Ai). We also 
note that in our numerical implementation, we chose a value of Ai proportional to h, i.e. Ai = Xh, where usually 
A = 1.0. 

We need to take special care if a particle leaves the computational domain (r^ > Tmax) or if it reaches the origin. In 
the first case we simply remove the particle from the integration scheme. When a particle reaches the origin, which 
operationally is signaled by r,; < 0, we "reflect" the particle by setting: 

n ^-n, (48) 

[Pr]^ ^ - [Pr], , (49) 

k -^ k (50) 



C. Interpolation and restriction 

In this section we explain how we calculate the stress-energy quantities on the finite-difference mesh from a given 
set of particles (restriction) , as well as how we interpolate the geometric quantities from the finite-difference mesh to 
the particles' positions. 

The values of the stress-energy quantities p, jV-, S are calculated on the mesh using equation ([2l|): 

N 

where fi are the single particle quantities. In our implementation, we use the specific kernel: 

TTz/ \ f l—ki'—'^il/^ '■ \fi—ri\<h ,-^s 

W{r,-^^) = [ '^ '\ : otherwise" (^^^ 

Similarly, to restrict the geometric quantities calculated on the mesh to the particles' positions we compute: 

F{n)=Y^F{r,)W{r,-n) (52) 

where, again, r^ is the position of the particle, rj are the grid points, and F is any of the coefficients which appear in the 
geodesic equations. These coefficients are generally products and quotients of metric functions and their derivatives. 
In order to calculate derivatives we use the standard 0{h^) finite-difference approximation: 

[F% = (F,+i - ^,_i) /(2/i) + 0{h^) (53) 

and then use equation ( [52|) to find an approximate value for F'{ri). 

D. Initial data 

To initialize the sets of particles which we evolve, we specify the particle distribution (number of particles per unit 
of areal coordinate) and the velocity distribution (specifically the number of particles per pr and /). This corresponds 
to a separable distribution function, f{r,pr,l): 

f{r,pr,l)^R{r)Pipr)Lil) (54) 

Moreover, instead of specifying P as a function of Pr we give P = P{pr) where Pr = Prja. This allows us to 
calculate the value of p* = ^Jm? +p^ + P /r"^ (assuming all the particles have the same rest mass w), and therefore 
[p]j) \jr]i and [S]i without a priori knowledge of the geometry. We can thus decouple the tasks of specifying initial 
conditions for the particles, and ensuring that the constraints are satisfied at i = 0. 

We use 1-dimensional Monte Carlo techniques applied to each of R{r), P{pr), L{1) to get a specific set of N particles. 
As mentioned above, the statistical error, in theory, should scale as 1/^/N. 

7 



IV. RESULTS 

All of the calculations discussed in this paper were performed with Nr — 257, N = 10^ and h — 0.078125. With 
this choice of parameters we ensured that the truncation error due to the finite-differencing of the field equations 
with mesh spacing h was of the same order of magnitude as the statistical error resulting from representation of the 
phase-space distribution with a finite number, N, of particles. We have observed Q that both types of error scale 
in the expected way: the truncation error scales as 0{h), with the number of particles per cell fixed; the statistical 
error scales as l/\/iV for fixed h, and as 0{l/h) for fixed N. Once the two errors are of comparable magnitude, in 
order to further decrease the overall error (truncation plus statistical error) as 0{h), we have to increase the number 
of particles as A^ « l//i^. This scaling behavior makes it very costly to substantially reduce the level of numerical 
error in the results presented here. 

Our critical solutions (solutions sitting at the threshold of black hole formation) were found by performing bisection 
searches using the total rest mass, Mq: 

N 

Mo ^^m^ = Nm (55) 

1=1 

as the tuning parameter. In particular, if M* is the critical parameter value, then configurations with Mq < M* 
(subcritical) will eventually disperse, while those with Mo > M* (supercritical) will form black holes. For any given 
critical search, we generally determined M* to a relative precision of about 4 x 10^^^. Tabic | provides a summary 
of the various families we have studied. 

We first focus on a specific family of initial conditions (Table fl Family (a)) and then summarize our observations 
for the remaining families. We thus consider an initial distribution defined by: 

R{r) = ^2g(-(r-r„)VA^) q^^^ (gg) 

P(p^) = e(-(p^-P-)VAp/) (57) 

L(0 = e(-('-'°)'/^'')e(0 (58) 

where, again, pr — Pr/cii and O is the step function. We take Tq — 5, A^ = 1, pro = 0, Ap^ — 2, Ig ^ 12, A/ = 2, 
and refer to the data as "almost time symmetric" (ATS) since the gaussian for the radial momentum is centered at 
Pr = 0. The critical parameter for this family is M* « 1.3 as shown in Fig 0. This figure also shows that the smallest 
black hole formed has finite mass and that the transition is therefore of Type I, in agreement with the observations 
of Rein et a/ |l| . 

In Fig. we show a few snapshots of the evolution of a{t,r) = dta{t,r) resulting from initial data which is close 
to criticality but which eventually disperses. At early times, a(t,r) oscillates, but for 100 ^ t < 200 it appears to 
approach 0. In the last snapshot (t — 234) we observe that d(t, r) has become negative, corresponding to dispersal of 
the particles. We also observe similar behavior for the time derivatives a and /3. 

In order to better see how small a{t,r) becomes, we show in Fig. Q details of a{t,r) at i = 156 for three different 
sets of particles sampled from the same initial distribution function. The difference between the solutions obtained 
with different sets give us an estimate of the statistical error, A5(a), in the calculation. From the figure we can see 
that for the most part, |As(d)| « |d|. This is not the case for r > 5.5 where the three calculations all seem to indicate 
a specific non-zero value for d; however, we suspect that the amplitude of this feature may decrease if we tune closer 
to the critical solution, and if we use greater resolution. More importantly, the region r ^ 5.5 accounts for only about 
5-10% of the mass of the near-critical configuration. 

These results are thus consistent with d — > in the critical regime, although more definitive proof would require 
significantly higher particle numbers, as well as higher finite-difference resolution. If we accept that the metric 
coefficients become independent of t in the critical regime, then the critical spacetime is stationary. If, in addition, 
the vector N"- — {d/dt}"' is orthogonal to the spatial hypersurfaces, then the spacetime is static, and the shift function, 
(3{t, r), must vanish. In Fig. |g we show the evolution of the shift function. During the period when the time derivatives 
of the metric coefficients are close to zero, the shift function /3(i, r) is also close to zero, in the sense that |A5(/3)| « |/3|. 
We thus have evidence that the critical solution in this case is characterized by a static geometry. 

We have also observed that in the critical regime the total current density, jr, tends to zero. This must be the 
case if the spacetime is static. However, as shown in Fig. |6|, the S^r component of the stress energy tensor is non- 
zero near criticality. This means that although on average there are the same number of particles with positive 
(outward-directed) and negative (inward-directed) Pr, the mean value of \pr\ does not vanish. 

As is typical of Type I critical solutions, as we tune Mq — > M*, the dynamical solution spends more and more time 
"close" to the putative static solution, and we expect to find power-law scaling of the time, r, (the "lifetime" of the 



near-critical configuration) spent in the critical regime as a function of In \Mo — M*\. Specifically, we expect the static 
critical solution (or solutions, since we are unable to demonstrate convincingly that the model has a unique critical 
solution, up to trivial rescalings) to possess exactly one unstable mode in perturbation theory, with an associated 
Lyapounov exponent which is simply the reciprocal of the scaling exponent, ct, in the lifetime scaling law: 

T r^-a\n\Mo~Mo'\ (59) 

Fig. shows a plot oi t = t — tc versus InjAfo — M*\ where t is the total time that the particles in the solution 
generated with parameter Mo are localized within r = Tq = 6, and tc is the same quantity for the solution closest 
to criticality. We show results from calculations using three different sets of particles and the same gaussian family 
previously discussed. Using the residual scaling freedom in the model (the equations of motion are invariant under 
t — > Ki, r — > Kr for arbitrary k > 0), we have also normalized each critical solution to have unit ADM mass: 

r^r/Ar(i*,r„,ax) (60) 

i^i/Af"(r,r„ax), (61) 

Here t* is defined to be the instant at which the time derivatives of the metric components are closest to zero for the 
solution closest to criticality. M'^(i*,rmax) is then the value of the mass aspect function at i = t*, r = rmax, again 
for the most nearly-critical solution. We can see that there is a rough linear relation between the lifetime r of the 
near-critical configurations and In \Mo — M*\: 

T --(5.2 ±0.2) In I A/To -M/l (62) 

where the quoted uncertainty is an estimate of the statistical error. 

Qualitatively, then, our results are similar to what has been observed in other instances of Type I critical collapse [|5| . 
We have also found results similar to those just presented by using 4 other gaussian families (Table |, Families (b)-(e)), 
each with pro = —4 and with varying lo's of 3, 5, 7 and 12; Tq = 5, A^ = 1, Ap^ — 2 and A/ = 2, as for Family (a)). 
Specifically, in each case we find that the critical geometry appears to be static. We note that for smaller values 
of lo, the mass in the critical solution gets increasingly concentrated near r = 0, making accurate evolution with a 
uniform finite-difference grid more difficult. Finally, we have studied a family with the following initial single-particle 
distributions (Table |, Family (f)): 

R{r) ex (l - tanh ((r - ro)/Arf) e(r) (63) 

P{pr) oc (1 - tanh{{pr-pro)/ApJ) (64) 

L{1) ex (l - tanh {{I ~ Q/^if) 0{l)- (65) 

Here we took Vo = 5, Ar = 1, pro — —4, Apr — 2, lo — 7 and A/ = 2. For this data we also find evidence that as 
AIo — > M* the geometry becomes static. 

In Fig. H we show profiles of 2M{t* ,r)/r for all of the different families considered, each separate profile being 
selected from the corresponding period of near-critical evolution. Again, since different initial conditions set different 
overall length scales for the problem, we have normalized the results using the rescaling given by equations (|6(]|-|6l|) . 
We see that, after normalization, the peak of 2M(i*, r)/r is roughly at the same radial location, r = r* = 2.3. We also 
find that the better resolved a solution is, the closer it conforms to the best resolved solution (Table ||, Family (a)). 
This provides some indication that there may be a universal critical solution in this model (up to trivial rescalings, 
r — > Kr, t — > nt), but again, we would need better finite-difference resolution and many more particles to verify 
this conjecture. This figure also shows that the maximum of 2M(t*,r)/r has a value of approximately 0.76. This 
immediately shows that the critical solution is not one of the clusters considered in p6[ , since there are no equilibrium 
Einstein clusters with maximum 2M{r)/r larger than 2/3. 

We have also estimated a defined by equation ( p9[ ) for the different families described above. Fig. S shows the 
lifetime scaling measured for the various initial data sets, where the quoted uncertainty in each value of a is the 
standard deviation of the slope computed from a least squares fit. The values that we have obtained for the scaling 
exponents are also collected in Table ||. 

Finally, we are also interested in investigating the dependence of the critical solutions on the distribution of angular 
momentum. In Fig. |l^ we show r^ S°'ait* ,r) for the different families of initial data we have studied. Here 



N 



r^-S^. = r^ Y.^{r-r.)—--l^_ (66) 



(see equation (p^), and t* is defined as previously. We note that r'^S'^ait^r), is a dimensionless quantity which 
measures the square of the angular momentum of the distribution of particles. As in Fig. y, we have again rescaled 
the radial coordinate (and time) so that the critical configuration has unit ADM mass. We see that there is no obvious 
agreement of the profiles calculated from different families of initial data; clearly more work needs to be done in order 
to clarify the effect of initial angular momentum distributions on critical evolution in this model. 

V. CONCLUSIONS 

We have studied critical behavior at the threshold of black hole formation for coUisionless matter with angular 
momentum and have corroborated the findings of Rein et al [nl that the black holes which form at threshold in this 
model are of finite mass (Type I behavior). Further, our results indicate that for families with non-zero angular 
momentum, the critical solution has a static geometry, with non-zero radial particle momenta. We have also found 
evidence for a lifetime scaling law which is to be expected for Type I critical solutions, and have some indications 
of universality. In order to produce more definitive results using our current approach, we would need to employ 
many more particles and better finite-difference resolution. Since the critical behavior in this model does not appear 
to generate structure on arbitrarily small scales, it seems unlikely that adaptive methods, such as those used in [||, 
would be of much help here. Thus, it may be that the development of a finite-difference code to solve the Vlasov 
equation directly in phase space would be the best route to more accurate results. Perhaps most importantly, this 
should provide a technique with better-understood, and better-controllable, convergence properties. 
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Family 


Form of Initial Data 


Set No. 


a 


(a) 
(a) 
(a) 


Gaussian, 
Gaussian, 
Gaussian, 


lo = 12 ATS 
lo = 12 ATS 
lo = 12 ATS 


1 
2 
3 


5.1 ±0.2 
5.3 ±0.2 

5.2 ±0.2 


(b) 


Gaussian, 


lo = 2. 


1 


5.7 ±0.2 


(c) 


Gaussian, 


lo = 5 


1 


5.5 ±0.2 


(d) 
(d) 


Gaussian, 
Gaussian, 


lo = 7 
lo = 7 


1 
2 


5.0 ±0.2 
5.0 ±0.2 


(e) 


Gaussian, 


lo = 12 


1 


4.9 ±0.2 


(f) 


Tanh, lo ~ 


7 


1 


5.9 ±0.2 



TABLE L Summary of critical searches described in the text. Listed are family label, form of initial data, set number (for 
families where multiple, independent, A^'-particle representations of the initial distribution function were used) and computed 
lifetime-scaling exponents, a. See the text and particularly equations (jSq-pq) and (BSJm), respectively, for detailed definitions 
of the "Gaussian" and "Tanh" initial data. Also note that ATS stands for "almost time symmetric", as discussed in the text. 
The quoted error for all values of a is the estimated statistical error for the ATS data (see equation (p2)). 
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FIG. 1. Illustration of the smoothing kernel, W{r — 
(uniform) finite-difference mesh points with A^ — tj+i 



W(r-r-) 




. Here, n is the position of a particle; 
-j = constant = h 



rj-i, rj, rj+i 



are the 
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M^ 
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FIG. 2. Black hole mass as function of total rest mass, Mo, for the "almost time-symmetric", gaussian family of initial data 
described in the text (Table M, Family (a)). We observe that the smallest black hole has a finite mass (Type I transition) at 
a critical parameter M^ ~ 1.3 (dashed line). Computationally, we have tuned M^ to a relative precision of about 4 x 10~^^, 
which is typical of the critical surveys discussed in this paper. The discrete jumps in the black hole masses for AIo > M* 
reflect the discrete nature of the finite-difference grid. We have made no attempt to "interpolate" the location of the black hole 
horizon in the finite-difference mesh; hence our mass estimates will always satisfy Mbh = kAr = kh, for some integer k. 
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FIG. 3. Evolution of a from a marginally subcritical calculation (Mo < M^) using Family (a). Note that at intermediate 
times Idl ~ 0. 
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FIG. 4. Plot of a(156, r) from three separate Family (a) calculations using distinct initial particle sets (A'^ = 10^). The 
scatter in the displayed datasets gives a rough indication of the level of statistical error As(d) in the computations. The plot 
shows that, at least for r < 5.5, where 93% of the mass of the putative static cluster is located, there is little or no correlation 
between the three sets. Thus, any non-zero value of a in the critical limit may be attributable to finite- A?^ statistical fluctuations. 
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FIG. 5. Evolution of l3 from a marginally subcritical calculation using Family (a). The shift function apparently vanishes 
during the same period of time as does d. Therefore, during this interval, we have evidence that A''" — (d/dt)'^ is orthogonal 
to the hypersurfaces, and, thus, that the geometry is static. 
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FIG. 6. Evolution of S^r{t,r) from a marginally subcritical calculation using Family (a). During the static regime, S^r{t,r) 
is bounded away from zero, showing that although pr is zero on average, \pr\ is not. 
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FIG. 7. Illustration of scaling law for the lifetime of near-critical configurations. The quoted uncertainty for each value of a 
is the standard deviation of the slope, which has been computed using a least squares fit. 
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FIG. 8. Comparison of near-critical solutions for different families of initial data (Table uj. We see evidence for a universal 
profile 2M{t* ,r)/r, where t* (different for each family) is the instant when the temporal derivatives of the metric components 
are minimized, and r has been rescaled for each family so that all critical solutions have unit ADM mass. The maximum value 
of 2M{r)/r is about 0.76 showing immediately that the critical solution cannot be one of the clusters considered in llq], since 
there are no equilibrium Einstein clusters with maximum 2M{r)/r larger than 2/3. (Moreover, in contrast to the configurations 
studied here, all particles in an Einstein cluster are in circular orbits.) 
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FIG. 9. Scaling behavior for different families of initial data. We observe near-critical lifetime-scaling behavior for all the 
families we have studied, as expected for Type I solutions (static or periodic solutions, with one unstable mode in perturbation 
theory). The quoted uncertainty in o is given by the standard deviation of the least-squares slope. The axes ranges vary 
somewhat from sub-plot to sub-plot; the values shown for Family (d) are representative. 
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FIG. 10. r"^ S°' a{t* , r) for the different families in the near-critical regime. As described in the text, this quantity is a 
dimensionless measure of the squared-angular-momentum of the distribution. In contrast to Fig. H, we see no particular 
evidence of a universal profile here. However, considerably more resolution (both in h and N) is needed to accurately assess 
the impact of angular momentum on critical collapse in this model. 
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